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Abstract 



Dilute Bose gas with attractive interactions is considered at zero tem- 
perature, when practically all atoms are in Bose-Einstein condensate. The 
problem is addressed aiming at answering the question: What is the optimal 
trap shape allowing for the condensation of the maximal number of atoms 
with negative scattering lengths? Simple and accurate analytical formulas 
are derived allowing for an easy analysis of the optimal trap shapes. These 
analytical formulas are the main result of the paper. 
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I. INTRODUCTION 



Physics of dilute Bose gases is a topic of high experimental and theoretical interest [1-5]. 
One of main parameters, characterizing the properties of these systems, is the intensity of 
atomic interactions defined by the value of the scattering length. The latter, by means of 
the Feshbach resonance techniques, can be varied in a wide range, including the change 
from positive to negative scattering lengths (see review articles [6,7] and references therein). 
Inverting the sign of the scattering length from positive to negative means the change of 
effective atomic interactions from repulsive to attractive. 

The properties of Bose systems with repulsive and attractive interactions are principally 
different. Thus, the homogeneous Bose gas with attractive interactions is known to be unsta- 
ble for any interaction intensity [8-10]. This happens because, under attractive interactions, 
the sound velocity of the homogeneous Bose gas becomes imaginary and, respectively, the 
isothermal compressibility negative and the Bogolubov spectrum complex [8-10]. 

The situation is different for trapped atoms, for which attractive atomic interactions can 
be stabilized by the zero-point kinetic energy due to the trapping potential, provided the 
number of atoms does not exceed a critical value iV c . Several experiments demonstrated 
the existence of the Bose-Einstein condensate of trapped atoms with negative scattering 
lengths. Thus, Bose-condensed 7 Li with a negative scattering length was trapped and studied 
in experiments [11-13]. Bose condensates of 85 Rb, with negative scattering lengths, were 
obtained by varying the latter from positive to negative by means of the Feshbach resonance 
techniques [14,15]. 

We may note that the Feshbach resonance techniques were also used for studying the 
instability in the dynamics of Bose-condensed gases, which developed owing to an abrupt 
change between the repulsive weak-coupling to strong-coupling regimes. For instance, in 
the experiment [16], a sample of 85 Rb atoms was evaporatively cooled and condensed into 
a state of atoms with a positive but a very small scattering length of about 80 Ob, where 
as is the Bohr radius. Then, by means of a fast variation of the magnetic field close to the 
Feshbach resonance, the scattering length was rapidly increased to 1900 as- The dynamic 
instability, observed in this rapid increase of the positive scattering length has a very different 
physical origin [16] and will not be discussed here. In the present paper, we shall consider 
only equilibrium Bose-condensed gases with attractive interactions, hence, with negative 
scattering lengths. 

Defining the critical number of Bose-condensed atoms, N c , which could be loaded into a 
trap, requires to accomplish heavy numerical calculations, since there are no small parame- 
ters that could facilitate the calculational procedure. The critical number N c was calculated 
by different numerical methods for spherically symmetric traps [17-24] and for anysotropic 
traps [25-27]. The general understanding for achieving the better collapse-avoiding prop- 
erties is as follows. Attractive atomic interactions are stabilized by the zero-point kinetic 
energy due to the trapping potential. In order to minimize the effect of attractive interac- 
tions for a given number of atoms, it is useful to reduce the condensate density, i.e., to use 
traps with lower confinement. The zero-point energy, corresponding to the weakest confin- 
ing trap direction, is the one relevant for determining the stability limits. When the radial 
trap frequency was held fixed, a cigar-shaped trap had a higher critical number mostly be- 
cause doing so lowered the condensate density from what it would have been if the trap 



2 



was spherical. Eventually, the lowered confinement in the weaker direction resulted in too 
little zero-point energy to stabilize the atoms. Thus, there should be an optimal critical 
number as one of the trapping frequencies would be relaxed. If the optimal trap shape were 
considered under a fixed condensate density, then the spherical shape would be optimum. 

Although the overall physical picture, resulting from numerical calculations [17-27], is 
understandable, there are no simple analytical relations allowing for a not too complicated 
investigation of the attractive Bose-Einstein condensate stability under the given parameters 
of an anisotropic trap. It is the aim of the present paper to advance a novel theoretical 
approach for the problem, using which we derive approximate analytical formulas connecting 
the critical number of atoms iV c with the trap frequencies. Though the derivation requires 
the usage of some elaborated mathematics, the final expressions are quite simple. Being 
simple, our formulas are at the same time rather accurate. Comparing them with the known 
numerical calculations [17-27], we find a good agreement, with a deviation within the range 
of the order 10%. No such analytical formulas have been obtained earlier. The derived 
analytical expressions can serve as a convenient tool for understanding and estimating the 
basic relations between the trap shape and the maximal number of trapped atoms with 
negative scattering lengths. 

Throughout the paper the system of units is employed, where the Planck constant, H — 1, 
and the Boltzmann constant, ks — 1, are set to unity. 



II. CONDENSATE WAVE FUNCTION 

We consider a dilute weakly interacting Bose gas, corresponding to the inequality 
p|a s | 3 <C 1, in which p is the average atomic density and a s is the s-wave scattering length. 
For repulsive interactions, a s > 0, while for attractive interactions, a s < 0. The dilute 
Bose gas in the low-temperature limit becomes almost completely condensed, forming Bose- 
Einstein condensate with the number of atoms N being approximately equal to the total 
number of atoms N. The condensate wave function, as is known [1-5], satisfies the Gross- 
Pitaevskii equation 



77(r) = e rj(r) , (1) 

where m is atomic mass, U(v) is a trapping potential, and the interaction strength is 

% = 4vr ± . (2) 
m 

Equation (1) has the form of the stationary nonlinear Schrodinger equation representing the 
eigenvalue problem, with the eigenfunction i](r) and eigenvalue e. The confining potential 
is usually described by the harmonic potential 

U(r) = y £ rl , (3) 

whose effective frequencies uo a define the anisotropy of the trap. The condensate wave 
function is normalized to the number of condensed atoms, 
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I \ V (r)\ 2 dr = N . (4) 

Under the considered conditions, this number is close to the total number of atoms, N ~ N. 

The standard traps usually have the geometry of a cylinder or are almost cylindrical, 
which we assume in what follows. Then there are two characteristic trap frequencies, the 
transverse, or radial, frequency 

UJj_=UJ X = UJy (5) 

and the longitudinal, or axial, frequency u z . The anisotropy parameter 

^ (6) 

defines the actual trap shape. When A < 1, the trap is called cigar-shaped, or pencil-shaped. 
For A = 1, the trap is spherical. And when A > 1, one terms the trap as disk-shaped, or 
pancake-shaped. Related to the characteristic frequencies u± and u z , there are two oscillator 
lengths 

L = , h = ■ (7) 

One more typical length is defined through the geometric- average frequency uq, 

1 1 /3 

l = - , UJo = (oJ± OJz) ■ (8) 

From Eqs. (5) to (8), it is straightforward to notice the relations between the characteristic 
frequencies, 

^ = A^V = $3 , (9) 
and between the oscillator lengths, 

= ^ = ■ (10) 



Employing the cylindrical coordinates r = {r±,(p,r z }, where r± = ^Jr^ + r^, it is conve- 
nient to introduce the dimensionless variables 



The equality 



- r± - r ' z m\ 

r= — , z = r- ( n ) 



[No 

V(r) = J-jr ip(r,^p,z) (12) 



defines the dimensionless wave function ip(r,{p,z) of the radial variable r G [0, oo), angle 
variable (p G [0,2n), and of the longitudinal variable z G (— oo,+oo). The normalization 
condition (4) transforms to 



4 



J \ip(r, </?, z) \ 2 rdr dip dz — 1 . (13) 
With the notation for the dimensionless energy 

E = — , (14) 

the Gross-Pitaevskii equation (1) takes the form of the eigenvalue problem 

H NLS i> = Ei>, (15) 
in which the nonlinear Schrodinger Hamiltonian is 

H NLS = -^ + \(r 2 + \ 2 z 2 )+g\iP\ 2 . (16) 
Here 

„ 2 d 2 Id Id 2 d 2 
V = 1 1 1 

and the notation for the effective coupling parameter 

g = 4n^N (17) 

is introduced. 

Equation (15), with the Hamiltonian (16), clearly shows that the energy E is a function 
of g and A. Real- valued solutions for the eigenvalue E = E(g,X) exist not for all parameters 
g and A. For some values of these parameters, the eigenvalue E can become complex. The 
appearance of the imaginary part in the energy E, means that the state with this energy 
level becomes unstable, having a finite lifetime that can be estimated as 

1 

r 



[in E\ 

The boundary on the plane g — A, where real-valued solutions for E disappear, defines the 
critical line <? C (A), which, through the relation 

</ c (A) = 4vr ^ N c , (18) 

gives the critical number of atoms iV c . For the number of atoms N < N c , the eigenproblem 
(15) possesses well-defined solutions, with real-valued energies. But, if No exceeds N c , the 
energy becomes complex, which means that the system is unstable and desintegrates. This 
critical number follows from Eq. (18) yielding 

^/M = ^(A) Al/2 /U 

4tt \a s J An \aj V ; 

Or, according to Eq. (10), one has 

N C = <f> A* . (20) 

To find the critical line <? C (A), it is necessary to solve Eq. (15) for varied parameters 
g and A. Directly solving this nonlinear equation, with an external potential, requires to 
invoke cumbersome numerical calculations. It would certainly be desirable to obtain an 
explicit analytical expression for the critical line <? C (A), which would permit us to accomplish 
a straightforward investigation of Eqs. (19) and (20). 
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III. OPTIMIZED PERTURBATION THEORY 



The behaviour of the energy E, as a function of g and A, can be found by using the 
optimized perturbation theory. This theory advanced in Ref. [28], is based on introducing 
into the calculational procedure, such as perturbation theory or iterative algorithm, control 
functions, which govern the convergence of the calculational scheme. Control functions can 
be incorporated either into the initial approximation or into the resulting series by changing 
the variables. Then, to guarantee the convergence of the sequence of approximations, the 
control functions are defined by means of an optimization procedure. This results in the 
sequence of optimized approximants. The specific features of the optimized perturbation 
theory is its high accuracy and the possibility of its usage for the problems with no small 
parameters. This theory is nowadays widely employed for various problems (see the review- 
type papers [29,30] and references therein). In particular, it has been successfully applied for 
calculating the interaction-induced shift of the condensation temperature of dilute Bose gas 
[31-36], yielding very accurate results, well agreeing with numerical Monte Carlo simulations 
[37-41]. 

To solve Eq. (15), we may start with the zero-order Hamiltonian 



V 2 

H = - — 

2m 2 



+ I ( M V + X 2 V 2 Z 2\ (21) 

in which u and v are control functions. The corresponding zero-order energy is 

E Q = (2ra + |m| + 1) u + (k + ^ Xv , (22) 

in which n — 0, 1, 2, ... is the radial quantum number, m — 0, ±1, ±2, ... is the azimuthal 
quantum number, and k — 0, 1, 2, ... is the axial quantum number. Employing the Rayleigh- 
Schrodinger perturbation theory, we get a sequence of the approximants 

E j = E j (g,\,u,v). (23) 

For brevity, we do not include explicitly the dependence on the quantum numbers. The 
control functions are defined in each order through optimization conditions, e.g., through 
the variation 

SE j = ^-6u + ^-6v = 0, (24) 
ou ov 

which gives the functions 

Uj = Uj{g, A) , Vj = vj{g, A) . (25) 
Substituting these back into Eq. (23), we obtain the optimized approximants 

Ej = Ej{g,\,Uj,Vj) , (26) 

with the control functions (25). 

To simplify the following expressions, it is useful to introduce an effective coupling 
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x = 2gV\ I , 



(27) 



where the integral 



(I^IM^I 2 ) 



uy Xv 



(28) 



is defined through the matrix element of the wave functions ip corresponding to the Hamil- 
tonian (21). Explicitly, the wave functions are 



x 



ip(r,ip,z) 



2n\ u 



m|+l 



(n + \m\)\ 



1/2 



r |m| exp 



u 



X 



tmtp 



1/4 



At; 

— I exp 



A 

— vz 
2 



Hi, 



Xv z) , 



/27T \/Wk\ \ vr 

where L™(-) is an associated Laguerre polynomial and Hk( ) is a Hermite polynomial. Then 
the integral (28) takes the form 



1 2 



(n+ \m\)\ 2 k k\ 



x 2\m\ e -2x 





4 f 




dx 1 


Jo 



e- 2t2 Ht(t) dt , 



which is a number not depending on the control functions u and v. In the first order, we 
have the energy 



E 1 = 



P 



u + 



1 



u 



1 



1 r- 



(29) 



in which the notation 

p = 2ra+H + l, g=(2A; + l)A (30) 
is used. The variation of Eq. (29) defines the control functions by the optimization equations 



1 \[v 
1- ^ + ^x = , 



1 ux 

1- ^7 + 



. 



These equations are to be supplemented by the boundary conditions 



lim u = lim v — 1 

s->0 s^O 



(31) 



(32) 



requiring that, when atomic interactions are switched off, the solutions be recovered corre- 
sponding to a harmonic oscillator. By their meaning, the control functions play the role of 
effective frequencies, because of which they have to be positive. From Eqs. (31) it follows 
that the range of definition for the solutions u and v depends on the sign of atomic interac- 
tions. Thus, for repulsive interactions, when g > and x > 0, these solutions are to be in 
the interval 



< u < 1 , 0<i;<l 



(x > 0). 



(33) 
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But for attractive interactions, with g < and x < 0, one has 



u > 1 , w > 1 (x < 0) . (34) 

Equation (29), together with the optimization conditions (31), gives the optimized approx- 
imant E = E\, which we shall analyze in what follows. This optimized approximant can be 
written in the form 



u 4 \ v 

in which u and v are the solutions to Eqs. (31). 

The set of Eqs. (31) and (35) defines the control functions u and v and the optimized ap- 
proximant E as functions of the effective coupling (27) and quantum numbers (30). Explicit 
dependence on these parameters can be obtained in the limits of weak and strong coupling, 
as is demonstrated in the Appendix A. 

Of special interest is the behaviour of the ground-state level, which, in the case of attrac- 
tive interactions, becomes unstable before the higher excited levels. For the ground state, 
one has n = m = k = 0, hence p — 1 and q — A. The integral (28) is / = (27r)~ 3 / 2 . The 
optimization equations (31) reduce to 

! : u x 

1 ; + V^x = 0, 1 ? + = , (36) 

U Z V Z Ay/V 



where the effective coupling (27) is 



2^ 

x = . . (37) 



(2tt) 3 /2 • 

The optimized approximant (35) becomes 

E=l + ±(v + ±). (38) 
u 4 V v J 

Expression (38) represents the ground-state energy with a very good accuracy for the whole 
range of the interaction strength, from zero to arbitrary large values. Thus, for weak cou- 
pling, Eq. (38) gives an asymptotically exact value 

^ l + ~2 + wk <9 ^ ±0) - <39) 

The accuracy slightly diminishes for increasing interactions. However, even in the limit of 
infinitely strong repulsive interactions, Eq. (38) yields the limit 

E ~ 0.547538 (g\) 2/5 (g -> oo) , (40) 

which is only 2% different from the Thomas-Fermi limit that is known to be asymptotically 
exact for g — > oo. For large attractive interactions, when g < 0, there is no simple asymptotic 
form for g — > — oo, since the system becomes unstable at a finite critical value g c (ty- The 
consideration of attractive interactions is more complicated and, to find explicitly the critical 
line <7 C (A), we need to involve some more techniques. 
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IV. ATTRACTIVE ATOMIC INTERACTIONS 



Now the consideration of the previous section will be continued with the specification 
for attractive interactions, when g < and x < 0. First of all, let us note that real-valued 
solutions for the ground-state energy (38), with the control functions from Eqs. (36), exist 
only for x > x c , such that 

\x c \ < 1 (0 < A < oo) . (41) 
The proof of this inequality is as follows. Equations (36) can be rearranged to the form 

£ " W^5) = ■ < 42 > 

As far as the right-hand side of Eq. (42) is positive, it should be that 

1 - \x\y/v > . 

From here, \x\ < 1/y/v. Since, according to Eq. (34), v > 1, if x < 0, then \x\ < 1. The 
latter results in Eq. (41). 

Because the critical value is smaller than one, \x c \ < 1, for any anisotropy parameter 
A > 0, it is reasonable to study, first, the situation at asymptotically small \x\ < \x c \ < 1. 
Making the replacement x = —\x\, we may analyze the behaviour at small \x\ for the control 
functions as well as for the energy. Excluding v in the optimization conditions (36), we get 
the equation 

x 2 (u 2 - l) V — A (u 2 — l) 4 + A x A u 8 = (43) 

for the control function u. Vice versa, excluding u, we have Eq. (42) for the control function 
v. Energy (38) can be represented in terms of one of these functions, for instance, 



u 4 



(u 2 - If x 2 u 4 
+ 



x 2 u 4 (u 2 - iy 



(44) 



Using Eqs. (42), (43), (38) or (44), we can find the expansions for all functions in powers of 
\x\. The expansions for the control functions are 



and 



The expansion for the energy is 



u~J2 u k\%\ k (45) 

k=0 



v ~^v k \x\ k . (46) 

k=0 



E~j2c k \x\ k . (47) 

k=0 
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In Eqs. (45) to (47), j = 0, 1,2, The first several coefficients for these expansions are 

given in the Appendix B. 

Expansions (45) to (47) have sense for asymptotically small \x\ <C 1. The value x c (X), 
should be smaller than one, in agreement with Eq. (41). However, \x c \ is not necessarily 
asymptotically small. This can be illustrated by the spherically symmetric case, when A = 1. 
Then u = v , and Eqs. (31) reduce to just one equation 

xu 5/2 + u 2 - 1 = (A = l). (48) 

The latter equation possesses real solutions only for x > x c , such that 

x c = = -0.534992 , (49) 

which is related to g c = —4.212958. Note that the critical value (49) follows exactly from 
Eq. (48). 

In this way, we need more general expressions for the sought functions, but not merely 
their asymptotic expansions (45) to (47). This can be achieved by resorting to a method 
of an effective summation of power series. A very accurate method of summation, whose 
mathematical foundation is based on the self-similar approximation theory [42-46], is the 
method of self-similar factor approximants [47-49] , which was shown to be more general and 
accurate than the method of Pade approximants, the latter being just a particular case on 
the class of rational functions. 

The self-similar factor approximant of the j-th order for a power series, say for expansion 
(45), has the form 

u*=u l[(l + A t \x\r , (50) 
i=i 

in which Nj = j/2, when j is even, and Nj = (j + l)/2, with A 1 set to one, if j is odd. 
The coefficients and the powers rii are defined by the accuracy-through-order procedure, 
that is, by expanding the approximant (50) in power series with respect to \x\ and equating 
in that expansion the same-order terms with those of series (45). In the same way the 
approximants v * for another control function are constructed. And, similarly, for the energy 
we obtain 

£* = c o n(l + C*Mr • (51) 
i=i 

The factor approximants (50) and (51) extrapolate the validity of the asymptotic series (45) 
and (47) to the whole region of finite \x\. As has been shown [47-49], such an extrapolation 
works very well for the whole region of the variable, from zero to infinity, being of especially 
good accuracy for the region between zero and the values of the variable of order one. It is 
exactly the latter region 0<|a;|<|x c |<l, which we here deal with. 

Having in hands the analytical expressions for the sought functions u*, v*, and E*, we 
are already close to our aim of finding an analytical form of the critical line x c (X) and, 
respectively, g c (X). For this purpose, it is sufficient to investigate the behaviour of the 
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functions (50) and (51), with varying |x| from zero to the point \x c \, where real positive 
solutions stop existing. The critical line x c (X) defines, according to relation (37), the critical 
line for the coupling 



9c(X) 



(27T) 



3/2 



2V\ 



x c (X) 



(52) 



The latter, by definition (18), gives the critical number of atoms N c as a function of the 
anisotropy parameter A. 

However, as relations (19) and (20) show, the critical number of atoms is not expressed 
solely through the anisotropy parameter A, but also involves one of the characteristic lengths 
(7) or (8). So, N c is a function of two parameters, assuming that the scattering length is 
fixed. For instance, if we choose the transverse length l± as one of the parameters and define 



a±(X) = \g c (X)\ , 



then the critical number of atoms is 



N c = a ± (X) 



(53) 



(54) 



But if the axial length l z is chosen, then, with the notation 

a z(X) = ^ \g c {X)\ , 



we have the critical number 



N c = a z (X) 



(55) 



(56) 



Finally, fixing the length l , and using the notation 



we get the critical number 



A 1 / 6 

«o(A) = — \g c (X)\ , 



N c = a (X) 



(57) 



(58) 



These formulas demonstrate that N c depends on two parameters. Therefore, maximazing 
one of the critical couplings (53), (55), or (57) for finding a maximal N c implies a conditional 
variation under one of the characteristic lengths being fixed. 

To find the conditional maxima of iV c , we need to define the critical line x c (X), where 
positive solutions for functions (50) or (51) stop existing. To estimate the accuracy of the 
obtained critical values x c , we constructed the factor approximants (50) and (51) up to 
the 17-th order and compared the values of x c for A = 1 with the exact value (49). The 
comparison shows good numerical convergence of the factor approximants, the sequence {u*} 
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converging a little faster than {E*}, because of which the values of x c obtained from the 
former sequence are more accurate. We shall take this into account in what follows, obtaining 
x c from the most accurate sequence. The error of x c , found from the factor approximants up 
to the 5-th order, is of the order of 10%. After the fifth-order approximant, the error quickly 
diminishes to about 1%. We compared in detail the behaviour of x c (X) and, respectively, of 
g c (X), found from the factor approximants up to the fifth order, in the whole range of the 
anisotropy parameter A. This behaviour of g c (X) translates into the properties of the related 
critical couplings (53), (55), and (57). Thus, for the transverse coupling (53), the limiting 
behaviour at very small or very large anisotropy parameters is 

ol(A) ~ 1.253 A 1/2 (A -> 0) , 

a±(X) ~ 0.627 X' 1/2 (A -> oo) . (59) 

For the axial coupling (55), we find the limits 

a z {\) ~ 1.253 A (A -> 0) , 

a z (\) ~ 0.627 (A — > oo) . (60) 

And for the critical coupling (57), we get 

a (A) ^ 1.253 A 2/3 (A -> 0) , 

a (A) ~ 0.627 A _1/3 (A -> oo) . (61) 

Choosing between the critical lines x c (X), obtained from the factor approximants of 
different orders, we keep in mind our main aim of deriving sufficiently simple analytical 
expressions that would be convenient to study. The accuracy of the second- and third-order 
factor approximants are very close to each other, with the errors of the order of 10%. The 
accuracy of the fourth-order approximant is slightly better, but the expression for x c (X) is 
much more cumbersome. Therefore, we limit ourselves by the second-order approximant, 
which gives the critical line 

X J\) = ^— . (62) 

v ; 1 + 2A K J 

Respectively, the critical line for the effective coupling (52) is 

g c (X) = -(2vr) 3 / 2 ^ . (63) 

The higher-order versions of the critical lines, with the comparison of their properties are 
given in the Appendix C. From Eq. (53), we obtain the radial critical coupling, 

the axial critical coupling (55), 
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and the average coupling (57), 



These simple formulas allow us to accomplish a straightforward search for an optimal 
trap shape. As is explained above, the result depends on the constraint under which the 
maximal critical number N c is defined. Thus, fixing the radial trap frequency u±, hence, the 
radial length l±, and varying the anisotropy parameter A, one has to maximize the radial 
critical coupling (64), which gives 

max ol(A) = 0.443 (A = 0.5) . 

A 

This tells us that, under the given conditions, the cigar-shaped trap, with A = 0.5, will be 
able to confine the maximal number of atoms (54). 

If one chooses to fix the axial trap frequency, u z , hence the axial length l z , then one 
should look for the maximum of the axial critical coupling (65), which yields 

max a z {\) = 0.627 (A — > oo) . 

A 

This means that, under the chosen conditions, the disk-shaped trap can contain the maximal 
number of atoms (56). 

And, if one prefers to fix the average trap frequency o;o, hence the length / , then the 
average critical coupling (66) is to be maximized, which leads to 

max a (A) = 0.418 (A = 1) . 

A 

Under these conditions, the spherical trap is preferable, if one wishes to confine the maximal 
number of atoms (58). 

These conclusions are in agreement with the general physical picture discussed in the 
Introduction. Let us emphasize that if the optimal trap shape is considered under a fixed 
condensate density at the trap center 

p(0)~A^) 3/ \ 

that is, under a fixed u>o, then the spherically symmetric trap is optimal, in agreement with 
the above analysis. Note also that fixing the density p(0), or frequency u Q , is equivalent to 
fixing the condensation temperature 



where £(3) is the Rieman zeta function. 



N 1 1/3 



C(3)J 
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For the purpose of estimating the critical number of atoms with attractive interactions 
for a given trap with the fixed trap frequencies, that is, with the given characteristic lengths, 
it is convenient to express the critical numbers (54), (56), and (58) directly through the trap 
parameters. Then, from any of these expressions, we obtain 

C V2 |a.|(Z2 + Z2 + Z,r (67) 
It is convenient to define the dimensionless quantity 

^Vldfnr (68) 

where 1$ = (l x lyl z ) 1 ^ ■ Equation (68) can be rewritten in terms of the trap frequencies. 
Introducing the notation for the average inverse frequency u inv through the equality 



1 



we find 




(69) 



Formulas (67) to (70) can serve as a useful tool for a fast and simple estimations of the 
critical number of atoms, with a negative scattering length, which could be loaded into a 
given trap. 



V. CONCLUSION 

We have considered a dilute gas of trapped Bose atoms at zero temperature, when prac- 
tically all atoms are in the ground state of Bose-Einstein condensate. Atoms are interacting 
through attractive forces, corresponding to negative scattering lengths. The number of pos- 
sible atoms in such a Bose-Einstein condensate is limited by a critical number N c depending 
on the trap parameters. The main result of the present paper is the derivation of simple 
analytical expressions for estimating this critical number N c . The investigation of the de- 
rived formulas clearly shows that the number of atoms N c cannot be expressed through 
the sole parameter of anisotropy. Therefore, in order to answer the question "What is the 
optimal trap shape for keeping the maximal number of atoms with attractive interactions" , 
it is necessary to specify under what conditions one is looking for this maximum. Thus, if 
the radial trap frequency u± is kept fixed and the anisotropy parameter A is varied, then 
the cigar-shaped trap with A = 0.5 is optimal. But if the axial trap frequency u z is fixed, 
then the disc-shaped trap is optimal. And when the geometric-average frequency ujq is fixed, 
with the anisotropy parameter being varied, then the spherically symmetric trap becomes 
optimal. The spherical trap is also optimal, when the condensate density is kept fixed. In 
general, for estimating the critical number of atoms N c , it is convenient to use formulas (67) 
to (70), which are represented in the form valid for arbitrary trap frequencies u> x , cu y , and 
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u z , hence, for arbitrary trap lengths l x , l y , and l z . It is easy to check that these estimates 
are in good agreement with experiment. 

For instance, in the experiments [12,50] with 7 Li, having the negative scattering length 
a s = — 27.3ob = —14. 5 A, where a# is the Bohr radius, the trap with the frequencies 
uj x = 2ir x 150.6 Hz, u y = 2n x 152.6 Hz, and uo z = 2n x 131.5 Hz was used. This, for the 
mass rriu = 11.6 x 10~ 24 g, translates into the characteristic lengths l x = 0.309 x 10 -3 cm, 
l y = 0.307 x 10~ 3 cm, and l z = 0.331 x 10~ 3 cm. In these experiments [12,50] the critical 
number was estimated being between 600 and 1300 atoms, that is, on average, N c ps 950. 
Our formula (67), for this case, gives iV c = 910. 

In the experiments [14,15] with 85 Rb, the scattering length could be varied by means 
of the Feshbach resonance techniques. The trap had the frequencies u x = 2tt x 17.24 Hz, 
uj y = 2n x 17.47 Hz, and uo z = 2tt x 6.80 Hz, that is uq = 2n x 12.7 Hz. The dimensionless 
parameter (70) was found to be equal to N c \a s \/l = 0.46 ± 0.07. From our formula (70) for 
this trap, we find N c \a s \/l = 0.47. 

In this way, the derived formulas allow one to quickly get rather accurate estimates for 
the critical number of atoms. The advantage of these formulas is their simple analytical 
representation. Therefore, to find iV c and to optimize the trap shape, one does not need to 
plunge into lengthy and complicated numerical calculations, but it is sufficient to employ 
the derived simple analytical formulas. 
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Appendix A. Properties of Optimized Approximants 

The control functions and the optimized approximant for the dimensionless energy, dis- 
cussed in Sec. Ill, have the following asymptotic behaviour in the weak-coupling and strong- 
coupling limits. 

In the weak-coupling limit, when g — > ±0, the effective coupling parameter (27) is such 
that x — > ±0. Then from the optimization equations (31), we find 

u ~ 1 - — x + P + 3q x 2 - 3p 2 + 16pg + 20g 2 ^ 
2p 8p 2 q 64p 3 q 2 

1 p + q 2 7p 2 + 20pq + 12q 2 „ 
v ~ 1 - — x + - — - x - — r . 

2g 4pg z 64jr<2" :> 

These expressions satisfy the boundary conditions (32). 

In the strong-coupling limit, when g — > oo and rr — > oo, we have the control functions 



V\ 1/5 x -2/ 5 + g 2 - 3p 2 ^-e/s + 3p 4 -pV-g 4 
v q J 5(pq 3 y/ 5 5pq 



7 2 \ 2/5 „,„ 2(r> 2 -2a 2 ) ( a\ 2/5 _ 6/5 f 6q A - 4p 2 q 2 - p 4 _ 2 



-(f) ^ ^&)' 



op 2 



The strong-coupling limit for attractive interactions, for which x — > — oo, is not defined, 
since the solutions to Eqs. (31) become complex. 

The weak-coupling limit for the optimized energy (35) can be written as an expansion 

E ~ a + aix + a 2 x 2 + a 3 x 3 (x — > ±0) 

in which the coefficients are 

a 1 p + 2q (p + 2q\ 2 

Respectively, the strong-coupling limit for the energy can be represented as 

E ~ 6 x 2 / 5 + 6i x- 2 / 5 + 6 2 x" 6 / 5 + 6 3 x- 2 , 
where the coefficients are 

5 / 2 vi/5 2p 2 + g 2 

3p 4 _ 2p 2 g 2 + 2g4 2j9 6 - p 4 q 2 - 2p 2 g 4 + 2g 6 

2 ~ 20(p 2 g) 3 /5 ' 3 ~ 20p 2 g ' 

These expansions are valid for arbitrary energy levels. The quantum numbers enter through 
notation (30). 
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Appendix B. Expansions for Attractive Interactions 



For atoms with attractive interactions, it is convenient to use the relation x = —\x\ and 
to express all quantities as functions of \x\. In expression (45) for the control function u, 
the coefficients are 

1 1 + 3A (1 + 2A)(3 + 10A) 
«o = l, «i=2> ^ = -83T' ^ = ^ , 

2 + 20A + 48A 2 + 35A 3 15 + 336A + 1400A 2 + 2048A 3 + 1008A 4 

Ua ~ 128A3 ' Ub ~ 4096A^ ' 

(1 + A) (35 + 221A + 409A 2 + 231A 3 ) 
M6 ~ 1024A 1 ' 

Expansion (46) for the control function v has the coefficients 

1 1 + A 7 + 20A + 12A 2 
Vo = 1 > Vl = T\> V2 = ^V> V3 = 64A^ ' 

5 + 32A + 48A 2 + 20A 3 39 + 616A + 1800A 2 + 1792A 3 + 560A 4 

V4 ~ 128A 4 ' V5 ~ 409QX 5 ' 

35 + 192A + 350A 2 + 256A 3 + 63A 4 

V « = 5T2A^ • 

And in expansion (47) for the optimized energy, the coefficients are 

_2 + A 1 1 + 2A (1 + 2A) 2 

C °~^ - ' Cl__ 2' ° 2 ''^OT 1 ° 3 ~ 64A 2- ' 

1 + 8A + 16A 2 + 10A 3 _ 3 + 56A + 200A 2 + 256A 3 + 112A 4 

° 4 ~ 256A^ ' ° 5 ~ 4096A 5 ' 

(1 + A) 2 (5 + 22A + 21A 2 ) 
° 6 ~ 1024A 4 ' 

These expansions are written here for the ground-state level, which looses stability before 
other states with higher energies for x — * x c and g — > g c , when attractive interactions become 
too strong. This is easy to understand keeping in mind that with increasing \x\ or \g\, that 
is, decreasing x and g, all energy levels move down. Consequently, the lowest energy level 
will be the first to touch zero and then become complex, provoking the system instability. 
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Appendix C. Properties of Critical Lines 



The critical line x c (X) is obtained as is explained in Sec. IV, from analysing the behaviour 
of the related self-similar factor approximants. On the line x c (X), the ground-state energy 
becomes complex- valued. Generally, for a fc-th order factor approximant, we get a k-th order 
critical line x[ k \\) and, respectively, the effective-coupling critical line 

9?\\) = ^x?\\). 

In the second order, we have 

x?\\) - 2X 



1 + 2A ' 

which is the critical line (62). In the third order, we get 

16A 



*W(A) = 



9 + 16A ' 

which is very close to xf\ And for the fourth order, we find 
.(4),^_ 8A(1 + 4A) 



^2(19 + 36A + 56A 2 + 64A 3 ) - 2 + 20A + 32A 2 ' 



The higher-order expressions x[ k \\) become quickly so much cumbersome that the whole 
idea of deriving simple analytical expressions, allowing for an easy analysis, looses grounds. 
This is why, in the main text of the paper the simplest forms for the critical lines (62) 
and (63), corresponding to the second-order approximations, were used. In addition, the 
functions g^ k \X) of all orders demonstrate rather similar behaviour. For instance, for the 
absolute values of the critical lines g( k \\), we obtain the following asymptotic behaviour in 
the limits of the small and large anisotropy parameter A. 

In the limit A — * 0, related to a very elongated cigar-shaped trap, we have 

|^ 2) (A)| ~ 15.7496 A 1/2 - 31.4992 A 3/2 + 62.9984 A 5/2 , 

\gi 3) (X)\ ~ 13.9997 A 1/2 - 24.8883 A 3/2 + 44.2458 A 5/2 , 
|#(A)| ~ 15.1278 A 1/2 - 33.3560 A 3/2 + 67.7766 A 5/2 . 
In the opposite limit A — > oo, which refers to a very plane disk-shaped trap, we obtain 

|^ 2) (A)| ~ 7.8748 A~ 1/2 - 3.9374 A~ 3/2 + 1.9687 A" 5/2 , 

|^ 3) (A)| ^ 7.8748 A- 1/2 - 4.4296 A~ 3/2 + 2.4916 A" 5/2 , 
|^ 4) (A)| ~ 7.8748 A~ 1/2 - 1.9687 \~ 3/2 + 1.7226 A" 5/2 . 
Note that the limits 

lim xi k) (X) = -1 , 



A^oo 



lim VAb?>(A)| = ^ 
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are the same in all orders. 

Looking for the optimal trap shape at a fixed condensate density, that is, under a fixed 
frequency u , we have to study the behaviour of the effective coupling 



4 fc) (A) 



A 1 / 6 

47T 



#(A) 



2tt 



4AV3 



x 



(k) 



(A) 



The maximum of this quantity defines the critical number (58). For all orders k, the max- 
imum of aQ k \\) happens at A m 1. Thus, A^ = 1, A® = 1.1, and A^ = 0.9. This means 
that, under a fixed condensate density, the spherical trap is optimal. 
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